# ============================================================
# Project Title: Running the hedonic treadmill: Lateral 
# insecurity and foreign real estate investment
# Replication Script: Experimental Results
# Date: 02/07/2026
# ============================================================

# ============================================================
# Clear environment and load packages
# ============================================================

rm(list = ls())

pacman::p_load(
  dplyr, haven, janitor, sjlabelled, jtools, interplot, ggplot2, gridExtra, grid,
  network, igraph, intergraph, tidyverse, readr, mediation, DiagrammeR, cobalt,
  psych, webshot2, gt, foreign,
  tidyr, showtext, scales, car, tibble
)

# ------------------------------------------------------------
# Paths / outputs
# ------------------------------------------------------------

WORKDIR <- "~/Desktop/Replication/"
OUTPUT_DIR <- file.path(WORKDIR, "replication_outputs")
if (!dir.exists(OUTPUT_DIR)) dir.create(OUTPUT_DIR)

setwd(WORKDIR)

# ------------------------------------------------------------
# Font setup
# ------------------------------------------------------------

# (You noted this works on your machine; keep as-is.)
font_add("Times New Roman", regular = "/System/Library/Fonts/Supplemental/Times New Roman.ttf")
showtext_auto()

# ============================================================
# Load data
# ============================================================

cces_data <- read_dta("unmatched.dta")

# Convert labelled variables to factors (optional, useful for consistency)
cces_data <- haven::as_factor(cces_data)

# Save converted version (replication artifact)
write.csv(cces_data, file.path(OUTPUT_DIR, "converted_from_dta.csv"), row.names = FALSE)

# ============================================================
# Clean and construct variables
# ============================================================

cces_data <- cces_data %>%
  mutate(
    
    # --------------------------------------------------------
    # Immigration attitudes
    # --------------------------------------------------------
    imm_att1 = case_when(CC24_323a == "Support" ~ 1, CC24_323a == "Oppose" ~ 0),
    imm_att2 = case_when(CC24_323b == "Support" ~ 0, CC24_323b == "Oppose" ~ 1),
    imm_att3 = case_when(CC24_323c == "Support" ~ 0, CC24_323c == "Oppose" ~ 1),
    imm_att4 = case_when(CC24_323d == "Support" ~ 1, CC24_323d == "Oppose" ~ 0),
    
    # (You used 1,3,4 — keep as-is)
    imm_att_index = rowMeans(across(c(imm_att1, imm_att3, imm_att4)), na.rm = TRUE),
    
    # --------------------------------------------------------
    # Economic perceptions
    # --------------------------------------------------------
    econ_pers = case_when(
      CC24_301 == "Gotten much better" ~ 2,
      CC24_301 == "Gotten somewhat better" ~ 1,
      CC24_301 == "Stayed about the same" ~ 0,
      CC24_301 == "Gotten somewhat worse" ~ -1,
      CC24_301 == "Gotten much worse" ~ -2,
      CC24_301 == "Not sure" ~ NA_real_
    ),
    
    # --------------------------------------------------------
    # Racial resentment
    # --------------------------------------------------------
    rr1 = dplyr::recode(CC24_441a,
                        "Strongly disagree" = 1,
                        "Somewhat disagree" = 2,
                        "Neither agree nor disagree" = 3,
                        "Somewhat agree" = 4,
                        "Strongly agree" = 5,
                        .default = NA_real_
    ),
    rr2 = dplyr::recode(CC24_441b,
                        "Strongly disagree" = 1,
                        "Somewhat disagree" = 2,
                        "Neither agree nor disagree" = 3,
                        "Somewhat agree" = 4,
                        "Strongly agree" = 5,
                        .default = NA_real_
    ),
    rr2 = 6 - rr2,  # reverse-code
    racial_resentment = rowMeans(across(c(rr1, rr2)), na.rm = TRUE),
    
    # --------------------------------------------------------
    # Demographics
    # --------------------------------------------------------
    age = 2022 - birthyr,
    age_group = case_when(
      age >= 18 & age <= 29 ~ 1,
      age >= 30 & age <= 39 ~ 2,
      age >= 40 & age <= 49 ~ 3,
      age >= 50 & age <= 59 ~ 4,
      age >= 60 ~ 5
    ),
    
    educ_cont = case_when(
      educ == "No HS" ~ 1,
      educ == "High school graduate" ~ 2,
      educ == "Some college" ~ 3,
      educ == "2-year" ~ 4,
      educ == "4-year" ~ 5,
      educ == "Post-grad" ~ 6
    ),
    
    gender_bin = case_when(
      gender4 == "Woman" ~ 1,
      gender4 == "Man" ~ 0
    ),
    
    # Income midpoints
    income_cont = case_when(
      faminc_new == "Less than $10,000" ~ 5000,
      faminc_new == "$10,000 - $19,999" ~ 15000,
      faminc_new == "$20,000 - $29,999" ~ 25000,
      faminc_new == "$30,000 - $39,999" ~ 35000,
      faminc_new == "$40,000 - $49,999" ~ 45000,
      faminc_new == "$50,000 - $59,999" ~ 55000,
      faminc_new == "$60,000 - $69,999" ~ 65000,
      faminc_new == "$70,000 - $79,999" ~ 75000,
      faminc_new == "$80,000 - $99,999" ~ 90000,
      faminc_new == "$100,000 - $119,999" ~ 110000,
      faminc_new == "$120,000 - $149,999" ~ 135000,
      faminc_new == "$150,000 - $199,999" ~ 175000,
      faminc_new == "$200,000 - $249,999" ~ 225000,
      faminc_new == "$250,000 - $349,999" ~ 300000,
      faminc_new == "$350,000 - $499,999" ~ 425000,
      faminc_new == "$500,000 or more" ~ 500000,
      TRUE ~ NA_real_
    ),
    
    # Perceived household income change
    income_change = case_when(
      CC24_302 == "Increased a lot" ~ 2,
      CC24_302 == "Increased somewhat" ~ 1,
      CC24_302 == "Stayed about the same" ~ 0,
      CC24_302 == "Decreased somewhat" ~ -1,
      CC24_302 == "Decreased a lot" ~ -2,
      TRUE ~ NA_real_
    ),
    
    # Perceived price change
    price_change = case_when(
      CC24_303 == "Decreased a lot" ~ -2,
      CC24_303 == "Decreased somewhat" ~ -1,
      CC24_303 == "Stayed about the same" ~ 0,
      CC24_303 == "Increased somewhat" ~ 1,
      CC24_303 == "Increased a lot" ~ 2,
      TRUE ~ NA_real_
    ),
    
    # Employment status (binary)
    employed = case_when(
      employ %in% c("Full-time", "Part-time") ~ 1,
      employ %in% c("Temporarily laid off", "Unemployed", "Retired", "Permanently disabled", "Homemaker", "Student") ~ 0,
      TRUE ~ NA_real_
    ),
    
    party_id = case_when(
      pid7 == "Strong Republican" ~ 1,
      pid7 == "Not very strong Republican" ~ 2,
      pid7 == "Lean Republican" ~ 3,
      pid7 == "Independent" ~ 4,
      pid7 == "Lean Democrat" ~ 5,
      pid7 == "Not very strong Democrat" ~ 6,
      pid7 == "Strong Democrat" ~ 7
    ),
    
    homeowner = ifelse(ownhome == "Own", 1, 0),
    white = ifelse(race == "White", 1, 0),
    
    urbancity = factor(as.character(urbancity)),
    
    # --------------------------------------------------------
    # Mediator
    # --------------------------------------------------------
    lat_ins = dplyr::recode(UCR420,
                            "Not difficult at all" = 1,
                            "A little difficult" = 2,
                            "Fairly difficult" = 3,
                            "Difficult" = 4,
                            "Extremely difficult" = 5,
                            .default = NA_real_
    ),
    
    # --------------------------------------------------------
    # Blame items (binary indicators)
    # --------------------------------------------------------
    imm_inflation = case_when(
      UCR415_9 == "selected" ~ 1,
      UCR415_9 == "not selected" ~ 0,
      TRUE ~ NA_real_
    ),
    
    foreign_inflation = case_when(
      UCR415_6 == "selected" ~ 1,
      UCR415_6 == "not selected" ~ 0,
      TRUE ~ NA_real_
    ),
    
    # --------------------------------------------------------
    # Dependent variable components
    # --------------------------------------------------------
    frei_att1 = dplyr::recode(UCR421,
                              "Strongly agree" = 3, "Agree" = 2, "Slightly agree" = 1,
                              "Neither agree nor disagree" = 0,
                              "Slightly disagree" = -1, "Disagree" = -2, "Strongly disagree" = -3,
                              .default = NA_real_
    ),
    frei_att2 = dplyr::recode(UCR422,
                              "Strongly agree" = 3, "Agree" = 2, "Slightly agree" = 1,
                              "Neither agree nor disagree" = 0,
                              "Slightly disagree" = -1, "Disagree" = -2, "Strongly disagree" = -3,
                              .default = NA_real_
    ),
    frei_att3 = dplyr::recode(UCR423,
                              "Strongly agree" = 3, "Agree" = 2, "Slightly agree" = 1,
                              "Neither agree nor disagree" = 0,
                              "Slightly disagree" = -1, "Disagree" = -2, "Strongly disagree" = -3,
                              .default = NA_real_
    ),
    
    policy_attitudes_index = rowMeans(across(c(frei_att1, frei_att2, frei_att3)), na.rm = TRUE),
    
    # --------------------------------------------------------
    # Treatment variable
    # --------------------------------------------------------
    treat = case_when(
      UCR419_text_split == "See UCR419 text" ~ 1,
      UCR419_text_split == "No UCR419 text" ~ 0
    )
  ) %>%
  filter(!is.na(econ_pers)) %>%
  mutate(inputstate = as.factor(inputstate)) %>%
  mutate(
    econ_insecurity_gap = price_change - income_change,
    econ_insecurity_bin = case_when(
      price_change > income_change ~ 1,
      TRUE ~ 0
    )
  )

# ============================================================
# Analysis dataset (complete cases)
# ============================================================

med_data <- cces_data %>%
  dplyr::select(
    treat, lat_ins, policy_attitudes_index,
    age_group, gender_bin, educ_cont, homeowner, income_cont,
    party_id, white, racial_resentment, imm_att_index,
    econ_pers, econ_insecurity_gap, inputstate, region, urbancity,
    imm_inflation, foreign_inflation, frei_att1, frei_att2, frei_att3,
    income_change, price_change, employed, econ_insecurity_bin
  ) %>%
  na.omit() %>%
  mutate(
    blame_score = rowSums(across(c(imm_inflation, foreign_inflation)), na.rm = TRUE),
    blame_attribution = if_else(imm_inflation == 1 | foreign_inflation == 1, 1, 0, missing = NA)
  )

# Optional: save analysis dataset used for replication
write.csv(med_data, file.path(OUTPUT_DIR, "analysis_dataset.csv"), row.names = FALSE)

# ============================================================
# 4. Construct validity and sample balance
# ============================================================

# Balance test

library(cobalt)

# Balance formula (same covariates you listed in your table)
bal_formula <- treat ~
  age_group + gender_bin + educ_cont + homeowner + income_cont +
  party_id + white + urbancity + employed + region

bal_obj <- bal.tab(
  bal_formula,
  data = med_data,
  estimand = "ATE",
  un = TRUE
)


# Extract standardized mean differences
balance_df <- bal_obj$Balance %>%
  as.data.frame() %>%
  tibble::rownames_to_column("Covariate") %>%
  dplyr::select(Covariate, Diff.Un) %>%
  dplyr::rename(SMD = Diff.Un) %>%
  dplyr::mutate(
    Type = dplyr::case_when(
      grepl("gender|white|homeowner|employed|region|urbancity", Covariate, ignore.case = TRUE) ~ "Binary",
      TRUE ~ "Continuous"
    )
  ) %>%
  dplyr::select(Covariate, Type, SMD)

print(balance_df)

# Correlations
cor.test(med_data$lat_ins, med_data$income_change, use = "complete.obs", method = "pearson")
cor.test(med_data$lat_ins, med_data$price_change, use = "complete.obs", method = "pearson")
cor.test(med_data$income_change, med_data$price_change, use = "complete.obs", method = "pearson")

# Cronbach's alpha
vars <- med_data %>% dplyr::select(lat_ins, income_change, price_change)
psych::alpha(vars)

# Factor analysis
vars_fa <- vars %>% na.omit()
fa_1 <- psych::fa(vars_fa, nfactors = 1, rotate = "oblimin", fm = "ml")
fa_2 <- psych::fa(vars_fa, nfactors = 2, rotate = "oblimin", fm = "ml")

print(fa_1$loadings, digits = 3, cutoff = 0.3)
print(fa_2$loadings, digits = 3, cutoff = 0.3)

# ============================================================
# 5. Main models — direct effects
# ============================================================

model_ate_latins_simple <- lm(lat_ins ~ treat, data = med_data)
summary(model_ate_latins_simple)

model_ate_frei_simple <- lm(policy_attitudes_index ~ treat, data = med_data)
summary(model_ate_frei_simple)

# Interaction models (treat x homeowner)
model_latins_interact <- lm(lat_ins ~ treat * homeowner, data = med_data)
summary(model_latins_interact)

model_frei_interact <- lm(policy_attitudes_index ~ treat * homeowner, data = med_data)
summary(model_frei_interact)

# ============================================================
# 6. Main models — mediation (state FE)
# ============================================================

med_model <- lm(lat_ins ~ treat + age_group + homeowner + white +
                  racial_resentment + party_id + gender_bin + educ_cont +
                  urbancity + income_cont + imm_att_index + blame_attribution +
                  price_change + income_change + employed + inputstate,
                data = med_data)

out_model <- lm(policy_attitudes_index ~ lat_ins + treat + age_group + homeowner + white +
                  racial_resentment + party_id + gender_bin + educ_cont +
                  urbancity + income_cont + imm_att_index + blame_attribution +
                  price_change + income_change + employed + inputstate,
                data = med_data)

summary(med_model)
summary(out_model)

set.seed(1234)
med_out <- mediation::mediate(
  model.m = med_model,
  model.y = out_model,
  treat = "treat",
  mediator = "lat_ins",
  sims = 2000
)
summary(med_out)

# ============================================================
# 7. Mediation — region FE (alternate spec)
# ============================================================

med_model2 <- lm(lat_ins ~ treat + age_group + homeowner + white +
                   racial_resentment + party_id + gender_bin + educ_cont +
                   urbancity + income_cont + imm_att_index + blame_attribution +
                   price_change + income_change + employed + region,
                 data = med_data)

out_model2 <- lm(policy_attitudes_index ~ lat_ins + treat + age_group + homeowner + white +
                   racial_resentment + party_id + gender_bin + educ_cont +
                   urbancity + income_cont + imm_att_index + blame_attribution +
                   price_change + income_change + employed + region,
                 data = med_data)

set.seed(1234)
med_out2 <- mediation::mediate(
  model.m = med_model2,
  model.y = out_model2,
  treat = "treat",
  mediator = "lat_ins",
  sims = 1000
)
summary(med_out2)

# ============================================================
# 8. Sensitivity analysis — full sample
# ============================================================

set.seed(1234)
med_out_for_sens <- mediation::mediate(
  model.m = med_model,
  model.y = out_model,
  treat = "treat",
  mediator = "lat_ins",
  sims = 1000
)

sens <- mediation::medsens(
  med_out_for_sens,
  rho.by = 0.1,
  effect = "indirect",
  sims = 1000
)

summary(sens)

pdf(file.path(OUTPUT_DIR, "fig_sens_main.pdf"), width = 6, height = 5)
plot(
  sens,
  sens.par = "R2",
  xlab = expression(R[Y]^2),
  ylab = expression(R[M]^2),
  main = expression("ACME(" * R[M]^2 * "," * R[Y]^2 * ")")
)
dev.off()

# ============================================================
# 9. Mediation by housing tenure (heterogeneous ACME)
# ============================================================

med_model1 <- lm(
  lat_ins ~ treat * homeowner + age_group + white + racial_resentment +
    party_id + gender_bin + educ_cont + urbancity + income_cont +
    imm_att_index + blame_attribution + price_change + income_change +
    employed + inputstate,
  data = med_data
)

out_model1 <- lm(
  policy_attitudes_index ~ lat_ins * homeowner + treat + age_group +
    white + racial_resentment + party_id + gender_bin + educ_cont +
    urbancity + income_cont + imm_att_index + blame_attribution +
    price_change + income_change + employed + inputstate,
  data = med_data
)

set.seed(1234)
med_out_homeowners <- mediation::mediate(
  model.m = med_model1,
  model.y = out_model1,
  treat = "treat",
  mediator = "lat_ins",
  covariates = list(homeowner = 1),
  sims = 2000
)

med_out_nonhomeowners <- mediation::mediate(
  model.m = med_model1,
  model.y = out_model1,
  treat = "treat",
  mediator = "lat_ins",
  covariates = list(homeowner = 0),
  sims = 2000
)

summary(med_out_homeowners)
summary(med_out_nonhomeowners)

# Helper to extract ACME row (kept your approach)
get_acme_row <- function(med_obj, label) {
  s <- summary(med_obj)
  est   <- as.numeric(s$d0)
  ci_lo <- as.numeric(s$d0.ci[1])
  ci_hi <- as.numeric(s$d0.ci[2])
  data.frame(group = label, estimate = est, lower = ci_lo, upper = ci_hi)
}

df_acme <- bind_rows(
  get_acme_row(med_out_homeowners, "Homeowners"),
  get_acme_row(med_out_nonhomeowners, "Non-homeowners")
)

df_sim <- bind_rows(
  data.frame(group = "Homeowners", effect = med_out_homeowners$d0.sims),
  data.frame(group = "Non-homeowners", effect = med_out_nonhomeowners$d0.sims)
)

# Plot settings: your preferred deep colors
old_blue <- "#2171B5"  # Homeowners in your note? (you mapped Non-homeowners -> blue earlier)
old_red  <- "#D9473C"

df_acme$group <- factor(df_acme$group, levels = c("Non-homeowners", "Homeowners"))

pad <- 0.05
ylim_min <- min(df_acme$lower) - pad
ylim_max <- max(df_acme$upper) + pad

p_acme <- ggplot(df_acme, aes(x = group, y = estimate, color = group)) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.5, color = "gray40") +
  geom_point(size = 3.2) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.12, linewidth = 0.9) +
  scale_color_manual(values = c(
    "Non-homeowners" = old_blue,
    "Homeowners"     = old_red
  )) +
  scale_y_continuous(expand = expansion(mult = c(0.02, 0.05))) +
  coord_cartesian(ylim = c(ylim_min, ylim_max)) +
  labs(x = NULL, y = "ACME (Indirect Effect)", color = NULL) +
  theme_minimal(base_size = 13) +
  theme(
    text = element_text(family = "Times New Roman"),
    legend.position = "none",
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    axis.text.x  = element_text(size = 14, margin = margin(t = 0)),
    axis.title.y = element_text(size = 14, margin = margin(r = 6)),
    axis.text.y  = element_text(size = 12),
    plot.margin = margin(4, 6, 2, 6)
  )

ggsave(
  file.path(OUTPUT_DIR, "homeowner_plot.pdf"),
  plot = p_acme,
  width = 4.6, height = 3.2, dpi = 300, device = cairo_pdf
)

# Sensitivity (non-homeowners)
sens_non <- mediation::medsens(
  med_out_nonhomeowners,
  rho.by = 0.1,
  effect = "indirect",
  sims = 1000
)

pdf(file.path(OUTPUT_DIR, "sens_non.pdf"), width = 6, height = 5)
plot(
  sens_non,
  sens.par = "R2",
  xlab = expression(R[Y]^2),
  ylab = expression(R[M]^2),
  main = expression("ACME(" * R[M]^2 * "," * R[Y]^2 * ")")
)
dev.off()

# ============================================================
# 10. Diagnostics (VIF)
# ============================================================

vif_med_df <- as.data.frame(car::vif(med_model)) %>%
  rownames_to_column("Variable") %>%
  mutate(
    VIF_Mediator = if ("GVIF^(1/(2*Df))" %in% colnames(.)) {
      `GVIF^(1/(2*Df))`
    } else if ("GVIF" %in% colnames(.)) {
      GVIF
    } else {
      NA_real_
    }
  ) %>%
  dplyr::select(Variable, VIF_Mediator)

vif_out_df <- as.data.frame(car::vif(out_model)) %>%
  rownames_to_column("Variable") %>%
  mutate(
    VIF_Outcome = if ("GVIF^(1/(2*Df))" %in% colnames(.)) {
      `GVIF^(1/(2*Df))`
    } else if ("GVIF" %in% colnames(.)) {
      GVIF
    } else {
      NA_real_
    }
  ) %>%
  dplyr::select(Variable, VIF_Outcome)

vif_table <- full_join(vif_med_df, vif_out_df, by = "Variable")
print(vif_table)

write.csv(vif_table, file.path(OUTPUT_DIR, "vif_table.csv"), row.names = FALSE)

# ============================================================
# 11. Figures (descriptives + direct effects)
# ============================================================

# Relabel variables for plotting
med_data_plot <- med_data %>%
  mutate(
    lat_ins_label = factor(lat_ins,
                           levels = 1:5,
                           labels = c("Not difficult at all", "A little difficult", "Fairly difficult", "Difficult", "Extremely difficult")
    ),
    treat_label = factor(treat, levels = c(0, 1), labels = c("Control", "Treatment"))
  )

# 11a. Lateral insecurity distribution by treatment
p_lat_ins <- med_data_plot %>%
  count(treat_label, lat_ins_label) %>%
  group_by(treat_label) %>%
  mutate(perc = n / sum(n)) %>%
  ggplot(aes(x = lat_ins_label, y = perc, fill = treat_label)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_fill_manual(values = c("gray50", "steelblue")) +
  labs(x = "Lateral Insecurity", y = "Percentage of Respondents", fill = "Condition") +
  theme_minimal(base_size = 14, base_family = "Times New Roman") +
  theme(
    plot.title = element_blank(),
    axis.text.x = element_text(angle = 30, hjust = 1),
    legend.position = "right"
  )

ggsave(
  filename = file.path(OUTPUT_DIR, "lat_ins_dist.pdf"),
  plot = p_lat_ins,
  width = 8, height = 5, dpi = 300, device = cairo_pdf
)

# 11b. FREI item distributions
frei_long <- med_data %>%
  dplyr::select(frei_att1, frei_att2, frei_att3) %>%
  pivot_longer(cols = everything(), names_to = "item", values_to = "response") %>%
  mutate(
    item = factor(item,
                  levels = c("frei_att1", "frei_att2", "frei_att3"),
                  labels = c("Broad Item", "Policy Specific #1", "Policy Specific #2")
    ),
    response_label = factor(response,
                            levels = -3:3,
                            labels = c(
                              "Strongly disagree", "Disagree", "Slightly disagree",
                              "Neither", "Slightly agree", "Agree", "Strongly agree"
                            )
    )
  )

response_colors <- c(
  "Strongly disagree" = "#E69F00",
  "Disagree"          = "#56B4E9",
  "Slightly disagree" = "#009E73",
  "Neither"           = "#F0E442",
  "Slightly agree"    = "#0072B2",
  "Agree"             = "#D55E00",
  "Strongly agree"    = "#003366"
)

p_frei <- frei_long %>%
  count(item, response_label) %>%
  group_by(item) %>%
  mutate(perc = n / sum(n)) %>%
  ggplot(aes(x = item, y = perc, fill = response_label)) +
  geom_col(position = position_dodge(width = 0.9), width = 0.75) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_fill_manual(values = response_colors, name = "Response") +
  labs(x = "FREI Policy Item", y = "Percentage of Respondents") +
  theme_minimal(base_size = 16, base_family = "Times New Roman") +
  theme(
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.text.x = element_text(angle = 20, hjust = 1),
    plot.title = element_blank(),
    legend.position = "right",
    legend.title = element_text(size = 12),
    legend.text  = element_text(size = 10)
  )

ggsave(
  filename = file.path(OUTPUT_DIR, "fig_frei.pdf"),
  plot = p_frei,
  width = 9, height = 5, dpi = 300, device = cairo_pdf
)

# 11c. Direct effects plot
get_treat_effect <- function(model, outcome_label) {
  cf <- summary(model)$coefficients
  
  est <- cf["treat", "Estimate"]
  se  <- cf["treat", "Std. Error"]
  df  <- model$df.residual
  
  tcrit <- qt(0.975, df = df)
  ci_lo <- est - tcrit * se
  ci_hi <- est + tcrit * se
  
  tibble::tibble(
    outcome  = outcome_label,
    estimate = est,
    ci_lower = ci_lo,
    ci_upper = ci_hi,
    sig      = !(ci_lo <= 0 & ci_hi >= 0)  # TRUE if CI excludes 0
  )
}

ate_data <- dplyr::bind_rows(
  get_treat_effect(model_ate_latins_simple, "Lateral Insecurity"),
  get_treat_effect(model_ate_frei_simple,   "FREI Policy Preferences")
)

effect_colors <- c("TRUE" = "#D9473C", "FALSE" = "#2171B5")

p_ate <- ggplot(ate_data, aes(x = outcome, y = estimate, color = sig)) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.1, linewidth = 1) +
  scale_color_manual(values = effect_colors) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
  labs(y = "Average Treatment Effect", x = NULL) +
  theme_minimal(base_size = 16, base_family = "Times New Roman") +
  theme(
    axis.title.y = element_text(margin = margin(r = 10)),
    axis.text.x  = element_text(size = 14),
    legend.position = "none"
  )

ggsave(
  filename = file.path(OUTPUT_DIR, "dir_eff.pdf"),
  plot = p_ate,
  width = 6, height = 4, dpi = 300, device = cairo_pdf
)

# ============================================================
# End
# ============================================================
